Crystalline Order On Riemannian Manifolds With 
Variable Gaussian Curvature And Boundary 



Luca Giomi* and Mark Bowick^ 
Department of Physics, Syracuse University, Syracuse, New York 13240-1130, USA 

We investigate the zero temperature structure of a crystalline monolayer constrained to lie on 
a two-dimensional Riemannian manifold with variable Gaussian curvature and boundary. A full 
analytical treatment is presented for the case of a paraboloid of revolution. Using the geometrical 
theory of topological defects in a continuum elastic background we find that the presence of a 
variable Gaussian curvature, combined with the additional constraint of a boundary, gives rise to a 
rich variety of phenomena beyond that known for spherical crystals. We also provide a numerical 
analysis of a system of classical particles interacting via a Coulomb potential on the surface of a 
paraboloid. 



I. INTRODUCTION 

Crystalline structures are ubiquitous in nature because 
ordered close-packed configurations frequently minimize 
the interaction energy between the component units of 
the system. In soft condensed matter systems, where 
curved or fluctuating geometries are energetically acces- 
sible, the interplay between order and geometry has been 
the subject of much attention in the last two decades [1- 
9]. The existence of some preferred geometry, as for in- 
stance that arising from the growth of a crystalline mono- 
layer on a rigid substrate, may influence the nature of the 
allowed order and, on the other hand, the formation of 
particular ordered structure may lead to a deformation 
in the geometry of the substrate (i.e. shape transition 
[10, 11]). Self assembled systems such as colloidosomes 
[4] or thin films of block copolymers [12] are realizations 
of such non-Euclidean "soft" crystals. Protein subunits, 
which comprise the shells of spherical viruses, also pro- 
vide an example in which the mechanical properties of 
the systems are affected by the formation of crystalline 
aggregates on surfaces equipped with a non-zero Gaus- 
sian curvature. 

Some progress in understanding crystalline arrange- 
ments of particles interacting on a curved surface was 
achieved by the authors of [2, 13] in the context of the 
geometric theory of topological defects. The original in- 
teracting particle problem is mapped to a system of inter- 
acting disclination defects in a continuum elastic curved 
background. The defect-defect interaction is universal 
with the particle microscopic potential determining two 
free parameters, the Young modulus of the elastic back- 
ground and the core energy of an elementary disclination 
[14]. 

In flat two-dimensional space particles almost always 
pack in triangular lattices, unless the interaction poten- 
tial is carefully tuned to select some other lattice topol- 
ogy. The addition of non-zero Gaussian curvature to the 



problem introduces frustration in the sense that perfect 
planar crystalline order is now incompatible with the cur- 
vature of the surface. The most clear-cut example of 
such geometrical frustration is the case of spherical crys- 
tals [3, 8]. The geometrical frustration in terms of lattice 
topology is revealed by the Euler theorem on the sphere 
which relates the number of vertices (V), faces (F) and 
edges (E) of any convex polyhedron : V — E + F = x, 
with % = 2 in the case of the 2-sphere. This classical 
result of topology can be rephrased in a form that is 
particularly useful to describe the presence of defects in 
the lowest energy configuration of a curved crystal by 
defining a topological charge as the departure from the 
ideal coordination number of a planar triangular lattice: 
qi = 6 — Ci, with Ci the coordination number of the z— th 
vertex. The Euler theorem states that the total disclina- 
tion charge of any spherical lattice must be equal to six 
times the Euler characteristic x : 



qi = 6 X = 12. 



(1) 
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Although identical particles with hard core repulsion in a 
plane typically pack into a triangular lattice with 6— fold 
coordination, any triangulation of the sphere is required 
to have twelve 5— fold disclinations (provided we restrict 
ourselves to the energetically-preferred charge q = ±1 
disclinations) as well as an arbitrary number of 6— fold 
vertices. These twelve extra disclinations are the conse- 
quence of the geometrical frustration associated the the 
curvature of the sphere. 

Topological defects can appear in spherical lattices 
in the form of isolated 5— fold disclinations placed at 
the vertices of an icosahedron, as frequently happens 
in the case of small viral capsids, or grouped into 
one-dimensional arrays of dislocations (tightly-bound 
(5, 7)— fold disclination pairs). These arrays, also known 
as grain boundary "scars" , appear on the sphere when 
the ratio R/a (R radius of the sphere, a Euclidean lat- 
tice spacing) exceeds a particular threshold value approx- 
imately equal to 5. 

In this article we present an analytical and numerical 
study of the defect structure of a two-dimensional crystal 
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constrained to lie on the surface of a paraboloid of revo- 
lution. The parabolic geometry introduces two novel and 
important features compared to the sphere: 1 ) a variable 
Gaussian curvature and 2) the presence of a boundary. 
Both these features must be treated properly for a thor- 
ough theoretical understanding. 

The paper is organized as follows. In Sec. II we briefly 
review the geometrical approach of [2, 3, 8] in which 
the basic degrees of freedom are the defects themselves 
rather than the interacting particles and we derive the 
zero-temperature energy of a paraboloidal crystal. This 
formalism has the advantage of reducing the number of 
degrees of freedom as well as being rather universal in 
the sense that it applies to a broad class of interacting 
potentials. In Sec. Ill we discuss the results obtained 
by the numerical minimization of the potential energy 
of a system of classical charged particles interacting via 
a Coulomb potential on the surface of a paraboloid in 
the light of the geometrical approach. Sec. IV will be 
devoted to conclusions. 



II. THE GEOMETRICAL APPROACH 



where is the total number of disclinations on the 
boundary and TVj is the number of disclinations in the 
interior of the paraboloid (in the following we will re- 
serve the letter N for the total number of defects and 
V for the total number of vertices of the lattice). The 
elastic free energy of the crystal may be expressed in the 
form [2, 3, 8, 13]: 



F = F el + F c + F 



(5) 



where Fq is the free energy of the defect-free monolayer 
and F c is the contribution to the free energy due to the 
core energy of disclinations and is proportional to the 
total number of defects. The elastic energy F e i associated 
with the defect interaction can be expressed in the form: 



F 
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d 2 x d 2 y G 2 l{x, y)p(x)p(y), 



(6) 



where Y is the Young modulus for the planar crystal. The 
quantity p{x) has the meaning of the effective topological 
charge density: 



N 



(7) 



A. The Elastic Free Energy 

Let P be the two-dimensional paraboloid of revolution 
l R 3 described in parametric form by: 



(2) 



where h is the height of the paraboloid and R the maxi- 
mum radius. In the following we will call K = 2h/R 2 the 
normal curvature of the paraboloid at the origin. The 
metric tensor gij (with determinant g) and the Gaussian 
curvature K are given respectively by: 




9ij = 



K(r) 







(1 + K 2 r 2 Y 



(3a) 



(3b) 



Such paraboloidal surfaces serve as a good testing ground 
for exploring the effects of both variable Gaussian curva- 
ture and the presence of a boundary on the nature of 
crystalline order. From the topological point of view a 
paraboloid of revolution is equivalent to a disk. The Eu- 
ler characteristic is thus x = 1. The total topological 
charge, taking into account the preferred coordination 
on the boundary and in the interior, is given in this case 
by: 

N b Ni 

Q = £(4 - Ci) + 5^(6 - Ci) =6, (4) 



where 8{x, Xi) is Dirac delta function on the manifold: 
5(x,Xi) = g~ 1 ^ 2 Y[i$( x ~ x i) an d G2h(x,y) represents 
the Green's function for the covariant biharmonic opera- 
tor on P. 

The calculation of the effective free energy (5) can be 
simplified if free boundary conditions are chosen. An 
application of the second Green identity to Eq.(6) leads 
straightforwardly to the form: 



(8) 



where x( x ) is the solution of the inhomogeneous bihar 
monic equation: 

A 2 x (x) = Yp(x) 



x g 



with boundary conditions 



X {x) = xedF 
u^xix) =0 xedP 



(9) 

(10a) 
(10b) 



in which V* is the usual contravariant derivative in the 
metric <?y and vi is the i— th component of the tangent 
vector v perpendicular to boundary. If the parametriza- 
tion (2) is chosen, the normal vector u is simply given by 
g r /\gr\, with g r = d r x the base vector associated with 
radial coordinate r. The solution of Eq.(9) will then be: 



X (x) = J d 2 yG L (x,y)T(y) 



(11) 



where Gi(a;,y) is the Green's function of the covariant 
Laplace operator on P with Dirichlet boundary condi- 
tions 



AGl{x, ■) = S(x, •) xeP 
G L (x,-)=0 xedl 



(12) 
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FIG. 1: The function T a (r) for different values of k in the 
range 1-2. 



the hclicoid or Schcrk surfaces. We will take this point 
up again in Sec. IV. 

Inserting Eq. (15) into Eq. (14), the function T(x) can 
be expressed in the simple form: 



^P- = f E 9iG L (x,Xi) - T s (\x\) + U(x), 



(17) 



i=i 



where the first term represents the bare contribution of 
the defects to the energy density and the second corre- 
sponds to the screening effect of the Gaussian curvature. 
Explicitly: 



r s (|x|) = io g 



q e 



Vl + K 2 !' 2 



1 + Vl + K 2 r 2 



(18) 



where r = \x\ and 



and T(x) = A%(cc) is the solution of the Poisson problem: 
AT(x) = Yp{x), (13) 
which can be expressed in the Green form: 
T{x) 



l + Vi 



Y 



= J d 2 yG L (x 1 y)p{y) + U(x) 1 



(14) 



where U(x) is an harmonic function on P that enforces 
the Neumann boundary conditions (10b). The calcula- 
tion of the Green's function Gl(x, y) satisfying equation 
(12) can be reduced to the more familiar planar problem 
by conformally mapping the paraboloid P onto the unit 
disk of the complex plane (see Appendix A for a detailed 
explanation) : 



G L {x,y) = ^log 



z(x) - z(y) 



1 - z{x)z{y) 



(15) 



where z(x) = ge 1 ^ , a point in the unit disk of the complex 
plane, is the image of a point on the paraboloid under the 
conformal mapping. The conformal distance g is related 
to r by: 



g(r) = A 



.Vl + K 2 ?' 2 



1 + Vl + n 2 r 2 1 



(16) 



with A a scale factor which ensures that g(R) = 1. As 
explained in detail in Appendix A, the Green's function 
Gi(a;,y), and hence the entire elastic free energy, de- 
pends only on the coefficients of the first fundamental 
form of the surface (i.e. the metric tensor gij). Thus 
the elastic energy associated with the defect interactions 
and thus the crystalline order is an intrinsic property of 
the manifold and so is invariant under local isometries. 
This observation, which might appear obvious in the case 
of isometric surfaces such the Euclidean plane and the 
cylinder, is quite remarkable when applied to more so- 
phisticated isometric manifolds such as the catenoid and 



exp (Vl + k 2 R 2 ) 



(19) 



is a normalization constant depending on boundary ra- 
dius R and the ratio n. 

Fig. 1 shows a plot of the screening function T s (r) 
for different values of k in the range 1 — 2. As expected 
the contribution due to Gaussian curvature is maximum 
at the origin of the paraboloid and drops to zero at the 
boundary. 

The calculation of the harmonic function U(x) requires 
a little more effort. If the crystal was defect-free (or 
populated by a perfectly isotropic distribution of defects) 
the function U (x) would be azimuthally symmetric and 
constant on the boundary. By the maximum principle of 
harmonic functions, U(x) would be then constant on the 
whole manifold and dependent only on k and the radius 
R: U(x) = U Ki r. Such a constant can be determined by 
integrating Ax{x) = T(x) and imposing the boundary 



U, 




K,R . 



2 3 
r 



FIG. 2: The quantity U K] r as a function of R for different 
values of k in the range 1-6. 



4 



condition (10b) [9], This gives: 



2tt f R 
U Ki r = —J dr^/gT s (r), 



where A is the area of the paraboloid: 
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2tt 

3^ 



(i + K 2 R 2 y -l 



(20) 



(21) 



As shown in Figure 2, the value of U Kt R quickly ap- 
proaches the linear regime as the size of the radius in- 
creases: 



kR 



(22) 



Then for a defect-free configuration, the contribution of 
the boundary to the energy density is a constant off- 
set that persist even for large radii. In the presence of 
disclinations, on the other hand, the function xi x ) is 
no longer expected to be azimuthally symmetric and the 
harmonic function U (x) will not be constant throughout 
the paraboloid. Introducing the harmonic kernel H (x, y) 
such that: 



U(x) = -J d 2 yH(x,y)p(y), 



(23) 



the determination of U(x) requires the calculation of the 
function H(x, y) associated with the Green's function of 
the weighted biharmonic operator arising from the con- 
formal map of the paraboloid onto the unit disk of the 
complex plane. This problem has been intensively stud- 
ied in the past few years by the mathematics community 
because it is connected with the extension of the maxi- 
mum principle for the weighted biharmonic operator of 
the form Auj" 1 A (sec Hcdenmalm et al. [15] for a good 
survey on this topic). In the case of radial weights, as 
arises from the conformal mapping of any surface of rev- 
olution, the function H(x, y) can be calculated explicitly 
(see Shimorin [16]). For the sake of completeness we re- 
port an exact expression of the harmonic kernel H(x,y) 
in Appendix B. The physical understanding of the solu- 
tion (17), however, doesn't require the complete solution. 
As shown from the numerical data of Sec. Ill, the dis- 
tribution of the defects along the boundary is predomi- 
nately symmetric, and thus we expect the constant factor 
(23) to be the leading contribution from the boundary 
even in the general case. 

For the ground state energy (zero temperature limit) 
that interests us we must minimize the energy (6) with 
respect to both the positions and the total number of 
defects. Since the energy density (17) depends on the 
difference between a curvature term and the defect den- 
sity we expect the defects to arrange themselves so to 
approximately match the Gaussian curvature. A com- 
plete screening of the Gaussian curvature would yield a 
crystal with zero elastic energy at zero temperature. On 
the other hand the core energy associated with a generic 
number N of disclinations will be linear in N and so it 





FIG. 3: Energy density T 2 (a;) for two different defect dis- 
tributions. The energy density corresponding to an isoiated 
+1— disclination at the origin is shown on the ieft. The defect 
distribution used for the figure on the right has been obtained 
from the numerical minimization of the Coulomb energy of a 
system of V = 50 charges (see Sec. III). 



will favor the fewest defects possible. The structure of 
the crystal at zero temperature will be governed there- 
fore by the competition between the energy cost of creat- 
ing a defect and the compensating screening effect of the 
Gaussian curvature. Once the distribution of the defects 
is known the elastic energy can be easily calculated by 
integrating the (18) on the whole paraboloid. In Fig. 3 
we show an example of the energy density T 2 (x) corre- 
sponding to two different arrangements of defects on a 
paraboloid with (k = 1.6 and R = 1). 



B. Large Core Energies: Pyramidal Lattices 

In the regime of large core energies F c 3> F e i , the cre- 
ation of defects is strongly penalized and the lattice nec- 
essarily has the minimum number of disclinations allowed 
by the topology of the paraboloidal substrate. From 
symmetry considerations we might expect the optimal 
distribution of defects to consist of b +1— disclinations 
arranged along the boundary at the base vertices of a 
6-gonal pyramid and a 6-fold apex (of topological charge 
go = 6 — b) at the origin. The homogeneous boundary 
conditions adopted require the first term in Eq. (17) to 
vanish when Xi £ d¥. In the minimal energy configura- 
tion then, the system has the freedom to tune the total 
number of defects along the boundary to minimize the 
clastic energy Eq. (8) for any given value of the ratio k. 
This behavior is exclusive to manifolds with boundary 
and doesn't have any counterpart in crystals on compact 
surfaces like the sphere and the torus. In the following we 
will see how this minimization leads to properties which 
we believe to hold, in the most general sense, on any 
surface with boundary. 

We will label a pyramidal configuration by Yj,, where 
b denotes the number of base +1— disclinations. The co- 
ordinates (r, </>) of the vertices are given by: 



(0,any), R, 



2-Kk 



(24) 



l<k<b ) 

Using the Euler theorem one can show that it is possible 
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TABLE I: Number of vertices V for four different Yb, n families 
in the range b £ [3, 6] (tetrahedron, square, pentagonal and 
hexagonal pyramid) and n £ [1, 10]. 

to construct infinite families of polyhedra with the sym- 
metry group Cbv from the pyramidal backbone Yf,. The 
number of vertices is given by: 

V = \bn{n + 1) + 1, (25) 

where n is a positive integer which represents the number 
of edges (not necessarily of the same length) of the poly- 
hedron which separates two neighboring disclinations. In 
the following we will refer to these polyhedra with the 
symbol Fig. 4 illustrates two Y^ n lattices for the 

cases 6 = 4 and n = 7 (with V = 113), and 6 = 5 
and n = 10 (V = 276). In Table I we report the num- 
ber of vertices for the four simplest Yj, jra polyhedra for 
n 6 [1, 10]. By a numerical minimization of the energy 
Eq.(8) one can establish that the Yj, are indeed equilib- 
rium configurations for b G [3,5], for some range of the 
parameters k and R. The cases of b = 5,6 are partic- 
ularly significant because they are characterized by an 
equal number of defects (N = 6) of the same topological 
charge (q = 1). The two configurations will be associated 
therefore with the same core energy F c and this intro- 
duces the possibility of a structural transition between 
Y5 and Yq governed by the curvature ratio n and the 
boundary radius R. For fixed R and small values of n the 
6— fold symmetric configuration Yq is the global minimum 
of the free energy Eq.(6). For k larger than some critical 
value k c (R), however, the Yq crystal becomes unstable 
with respect to the 5— fold symmetric configuration Y5. 
A numerical calculation of the intersection point between 
the elastic energies of Y5 and Y& for different values of k 
and R allow us to construct the phase diagram shown in 
Fig. 5. The word "phase" in this context refers to the 
symmetry of the ground state configuration as a function 
of the geometrical system parameters k and R. 

The scenario depicted in Fig. 5 can be understood 
heuristically by imagining a system of spherically sym- 
metric equally sized subunits initially arranged on the 
surface of a planar disk (k = 0). The most efficient pack- 
ing of this system is clearly the one in which the subunits 
are arranged in a triangular lattice with six 3— fold sites 
on the boundary at the vertices of a hexagon. If now we 




FIG. 4: Two examples of Yb, n triangulations of the paraboloid 
(F4.7 on the top and Ys,io on the bottom). Plaquettes with 
disclinations are highlighted in red, for +1— disclinations, and 
green for +2— disclinations. 



slightly deform the disk into a low-curvature paraboloid 
(k > 0) we might expect the hexagonal configuration 
to persist for small values of k. When the deformation 
is more pronounced, however, the curvature at the origin 
will be enough to support the existence of a 5— fold vertex 
and the system will undergo a structural transition from 
the Yq to the Y5 phase. In principle, if we keep increasing 
the curvature we might expect the crystal to undergo a 
further transition to the Y4 phase. In this case, however, 
the core energy will also increase by a factor 4/3 and 
so this is not generally possible in the regime in which 
F c ^> F e i. For intermediate regimes (i.e. F c ~ F e i), 
Y5 — > Y4 and Y4 — * Y3 transitions are also possible. The 
critical value of the parameters k and R, however, is not 
universal and will depend in detail on the values of the 
core energy and the Young modulus. 
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FIG. 5: Phase diagram in the large core energy regime. For 
small k the lattice preserves the 6— fold rotational symmetry 
of the flat case. As the curvature at the origin increases the 
system undergoes a transition to the Y5 phase. 



C. Small Core Energies: Scars And Coexistence 

When the core energy F c is small, the elastic energy 
Eq.(6) can be lowered by creating additional defects. In 
this section we present a simple argument for the ap- 
proximate phase diagram of a paraboloidal crystal in this 
regime. A more detailed analysis of the minimal energy 
configurations following from the solution Eq.(17) will be 
reported elsewhere. 

Let us assume that a fivefold disclination is sitting at 
the point Xo = (ro, <f>o) of P. We can introduce a notion 
of distance on the paraboloid by setting up a system of 
geodesic polar coordinates (s, p) with origin at x . We 
expect that the stress introduced by the defect is con- 
trolled by an effective disclination charge inside a circular 
domain Cl of geodesic radius L: 



leff 



2ji 



dp / dsy/g K(s,cp), 



(26) 



where q = 7r/3 is the charge of the isolated defect and the 
integral measures the screening due to the total Gaussian 
curvature within the domain. The metric tensor and the 
Gaussian curvature of a generic Riemannian manifold can 
be expressed in geodesic polar coordinates in the form 
(see for example Do Carmo [17]): 



9ij 



K(s, ip) 



1 

G 



d 2 s VG 
G 



(27a) 



(27b) 



metric around the origin (0, p) yields: 

VG = s-lK s 3 + o(s 5 ). 
For small distance from the origin, Eq.(26) becomes: 

q eff =q+ dip dsd 2 s VG (28) 
Jo Jo 

= q- nK L 2 + o(L 4 ). (29) 

The right hand side of Eq. (29) is a very general expression 
for the effective disclination charge at small distance and 
doesn't depend on the embedding manifold. If a grain 
boundary is radiating form the original disclination, we 
expect the spacing between consecutive dislocations to 
scale like l/g e //, with a the lattice spacing [2]. When 
q e ft — > + the dislocation spacing diverges and the grain 
boundary terminates. Since the Gaussian curvature is 
not constant on P, the choice of the origin (i.e. the posi- 
tion of the central disclination along the grain boundary) 
affects the evaluation of q e f / . One can identify upper and 
lower bounds by observing that: 



maxK(r) = K(0) = k 2 

r 

min/v(r) = K(R) 



(1 + K 2 R 



2u2\2 



(30a) 
(30b) 



Unlike the case of surfaces with constant Gaussian curva- 
ture, we expect the phase diagram for paraboloidal crys- 
tals to consist of three regions separated by the curves: 



KM 



K K 



(31) 



where G = g v ■ g v . Furthermore, an expansion of the 



When L — L(K m i n ) — > + the effective disclination charge 
goes to zero and the distance between two consecutive 
dislocations diverges at any point on P. On the other 
hand if L — L(K max ) — > 0~, the disclination charge will 
prefer to be delocalized in the form of grain boundary 
scars. For L(K max ) < L < L(K m i n ) the paraboloid will 
be equipped with both regions where the Gaussian cur- 
vature is high enough to support the existence of isolated 
disclinations and regions where the screening due to the 
curvature is no longer sufficient and the proliferation of 
grain boundary scars is energetically favored. This leads 
to a three region phase diagram in which the regime of 
isolated disclinations is separated from the delocalized 
regime of scars by a novel phase in which both isolated 
disclinations and scars coexist in different parts of the 
paraboloid according to the magnitude of the Gaussian 
curvature. 

To compare this result with the structure of 
paraboloidal crystals obtained by numerical minimiza- 
tion of the Coulomb energy (see Sec. Ill) we need to 
measure the distance L c in terms of the lattice spacing 
a and rephrase Eq.(31) as a condition on a (or cquiva- 
lently on the number of vertices V). To do this we note 
that in order for the domain Cl to completely screen 
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FIG. 6: Defect phase diagram for a paraboloidal crystal of 
radius R = 1. The two critical lines that separate the iso- 
lated disclinations (ID) regime from the coexistence regime 
and the coexistence regime from the scar phase correspond 
to the solutions of Eq.(32) for Ko = K m ia and Ko — i^max, 
respectively. 



the topological charge of the shortest scar possible (i.e. 
5 — 7 — 5), the geodesic radius L has to be large enough 
to enclose the entire length of the scar. Calling I the 
geodesic distance associated with a single lattice spac- 
ing a, we will then approximate L ~ 3£. In the spirit 
of the local expansion Eq.(29), we can then approximate 
the neighborhood of the central point Xq with a spherical 
cap of radius 1 / y/Ko- With this choice we obtain: 



2 



1 



(32) 



which, combined with Eq.(31), gives an equation in a and 
k for fixed R: 



2 



1 



cos 



3\/3 



(33) 



The lattice spacing a can be approximately expressed as 
a function of the number of vertices of the crystal by 
dividing the area A of the paraboloid by the area of a 
hexagonal Voronoi cell of radius a/2: 



.4 



2 ' 



(34) 



The phase diagram arising from the solution of Eq.(32), 
together with the condition Eq.(34), is sketched in Fig. 
6 for the case R = 1. The two critical lines that sepa- 
rate the isolated defects (ID in the plot) regime from the 
coexistence regime and this one form the grain bound- 
aries phase correspond respectively to the solutions for 



K min and K 



The simplicity of the 



criteria used to derive Eqs.(32) and (34) doesn't allow 
us to predict the regions surrounding the critical lines 
with high numerical accuracy, but does provide a semi- 
quantitative picture of the novel phenomenology of de- 
fects in non-Euclidean crystals that is generally sup- 
ported by the numerical results presented in Sec. III. 



III. NUMERICAL EXPERIMENTS 

A. Energy Minima And Complexity 

In the following section we report the results of a nu- 
merical minimization of a system of V classical particles 
interacting via a Coulomb potential E = J2i<j ^/\ x i~ x j I 
on the surface of a paraboloid. The equilibrium config- 
uration arising from this optimization problem can be 
viewed as a direct realization of a paraboloidal crystal 
and thus provides a testing ground for our analytical re- 
sults. 

The determination of the equilibrium properties of 
complex systems is complicated by the rich topography 
of the energy landscape, with its many, often deep, lo- 
cal minima (valleys) separated by high barriers (passes). 
The number of local minima grows rapidly with system 
size, making it increasingly difficult, or impossible, to 
find the global minimum. 

The effort in solving a given global optimization prob- 
lem is described by computational complexity theory. 
Locating the global minimum for a potential energy sur- 
face belongs to the class of problems known as NP-hard, 
for which there is no known algorithm that is certain to 
solve the problem within a time that scales as a power of 
the system size. 

The Thomson problem [18-25] of finding the optimal 
configuration of V interacting charges on a 2-sphere rep- 
resents, in this context, a celebrated example of a hard 
optimization problem. The existence of novel arrays of 
topological defects in minimal energy configurations pro- 
vides further insight into the structure of the energy land- 
scape. Computer experiments on the Thomson problem 
indicate that, in the range 70 < V < 112, the number 
of local minima for each value of V grows exponentially: 
TV ~ 0.382 exp(0.04971/) [21]. This trend is believed to 
continue for larger values of V, making the determina- 
tion of the global minimum a formidable computational 
challenge. In the case of the paraboloid we believe the 
prefactor in this scaling law is larger due to the additional 
constraint of the boundary. 



B. Paraboloidal Coulomb Crystals 

To construct equilibrium lattice configurations we 
adopted a parallel implementation of Differential Evolu- 
tion (DE) [26] (see App. C for an introductory descrip- 
tion of the algorithm together with the parallelization 



8 



V 


V 


V 1 


\ n 
1 u 




Kin pv p*v 


i n 
1U 


U 


U 


4 


D 


oo.y44oooyoy / 4U10 


on 
zU 


U 


i 

4 


<-:• 
O 




1 7fi COni A 0Q l J77On7 

i /y.ozyi4ooo / /zy/ 


on 
oU 


U 


D 


Iz 


Iz 


/ion H/in7yi7Qt;Qn/in7 
4oy.U4y f 4 f o0oU4U / 


A A 

4U 


z 


- 


1 s 
lo 


10 


sis SQnn^o^cn/in^G 
olo.ooUUDzooU4UDy 


tin 


4 


1 

4 


Z4 


1 Q 

lo 


1 QOI C7QGQ/1 Q070 
lozl.o f ooy404oz f Z 


OU 


z 


i n 
1U 


OQ 

zo 


on 
ZU 


lO/IO OQnOOl /inQ7CQ 

iy4y.zoUzyi4Uo / oo 


/U 


o 
z 


1 c 


ZD 


o e, 
ZD 


07m ncntifinc: /i i oo i 

z /ui.yoyoouo4izzi 


80 
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17 


31 


29 


3581.110585181344 


90 


2 


16 


46 


26 


4588.364706108566 


100 


3 


15 


55 


27 


5722.503370970009 


150 


1 


30 


81 


38 


13323.70617345018 


200 


3 


35 


115 


47 


24173.215803305491 



TABLE II: Numerical data for twelve selected lattices. The 
quantities V q represent the number of vertices in the crystal 
with topological charge q. 



strategy). The initial pool of candidate solutions is gen- 
erated at the beginning of the simulation by randomly 
creating NP — 2QV configurations uniformly distributed 
over the whole search space {re [0, R]} ® {4> G [0, 2ir]}. 
The population is then evolved by 3 • 10 5 DE iterations 
on ten processors working in parallel. 

In Fig. 7 we show the Voronoi lattice and the Delau- 
nay triangulations for five selected systems up to V = 200 
particles. The lowest energy seen, together with the num- 
ber of n— fold vertices, for each one of these lattices is re- 
ported in Table III A. In all the systems observed discli- 
nations always appear clustered in cither grain bound- 
ary scars or dislocations with the exception of isolated 
+1— disclinations which appearing in the bulk as ex- 
pected from the curvature screening argument discussed 
in Sec. II. The complex aggregation of defects along 
the boundary together with the presence of negatively 
charged clusters indicates that the effect of the boundary, 
in the case of relatively small systems like the ones simu- 
lated, is more drastic than predicted by the homogeneous 
boundary conditions Eq.(lOa) and Eq.(lOb). Even in the 
computationally expensive case of V — 100, the distance 
between the origin and the boundary of the paraboloid 
is only four lattice spacings. In this situation we expect 
the distribution of particles along the boundary to play 
a major role in driving the order in the bulk. 

For larger systems, such as V = 200 (top of Fig. 7), 
the behavior of the particles in the bulk is less affected 
by the boundary and the crystalline order reflects more 
closely the free-boundary problem discussed in Sec. II. 
A comparison of the lattices in the first two rows of Fig. 
7, in particular, reveals substantial agreement with the 
scenario described in Sec. II B. For k = 0.8 and V = 200, 
the defects are all localized along the boundary with the 
exception of one length-3 scar in the bulk at distance 
r ~ 0.63 from the center. For k = 1.6, the pattern of 
defects in the bulk is characterized by the coexistence of 




FIG. 7: Voronoi lattices and Delaunay triangulations for five 
selected systems from numerical simulations with R = 1. The 
first row corresponds to V — 200 and k = 1.6, while the 
second row is for V = 200 and k = 0.8 (see Sec. IIIB for 
a discussion). In the bottom two rows V — 100, 80, 50 with 
k = 1.6. 



an isolated +1— disclination at the origin and a length- 
5 W— shaped scar displaced along a parallel one lattice 
spacing away from the central disclination. Apart from 
the evident difficulty in comparing the structures of small 
systems with those predicted from continuum elasticity 
theory, this behavior is consistent with the simple pic- 
ture sketched in the phase-diagram of Fig. 6. The local 
5— fold symmetry at the origin of the k = 1.6 config- 
uration, compared with 6— fold symmetry for k = 0.8, 
suggests, as in the case of spherical crystals [27], that 
the complicated structure of defect clusters appearing in 
large systems is the result of the instability of the sim- 
pler Yb,n configurations from which they partially inherit 
their overall symmetry. A more accurate numerical veri- 
fication of our theory remains a challenge for the future. 
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The symmetry of the configurations presented in Fig. 
7 deserves special attention. As for any surface of revolu- 
tion, the circular paraboloid belongs possesses the sym- 
metry group 0(2) of all rotations about a fixed point 
and reflections in any axis through that fixed point. Any 
given triangulation of the paraboloid may destroy the full 
rotational symmetry completely or just partially, leaving 
the system in one of the following two subgroups: the 
pyramidal group C nv or the reflection symmetry group 
C s . In general we found the latter symmetry group for 
system sizes up to V = 200 particles. The symmetry for 
larger system sizes is under investigation. 



IV. CONCLUSIONS AND DISCUSSIONS 

We have analyzed both analytically and numerically 
the structure of a two-dimensional paraboloidal crystal as 
a specific realization of the class of crystalline structures 
on two-dimensional Ricmannian manifolds with variable 
Gaussian curvature and boundary. Using the geometrical 
approach developed in [2] we found that the presence of 
variable Gaussian curvature, combined with the presence 
of a boundary, gives rise to a rich variety of phenomena 
which we believe to be generic. The freedom of tuning the 
optimal number of disclinations to place on the boundary, 
where the presence of defects doesn't cost additional en- 
ergy, allows the crystal to undergo a structural transition 
controlled by the value of the maximum Gaussian curva- 
ture of the surface, in the regime where the creation of 
defects is energetically disfavored and disclinations are 
isolated in the crystal and far apart from each other. 
When the core energy of the defects is much smaller than 
the elastic energy and disclinations are allowed to prolif- 
erate, the presence of a variable Gaussian curvature is 
responsible for the existence of an intermediate phase in 
which both isolated defects and grain boundary scars co- 
exist in the crystal in different regions according to the 
magnitude of the curvature. 

By conformally mapping the paraboloid into the unit 
disk of the complex plane, we showed furthermore that 
the elastic energy of the crystal only depends on the co- 
efficients of the first fundamental form of the surface and 
thus is an intrinsic property of the manifold which is 
invariant for local isometrics. This property, which is 
somehow already contained in the linearity of the elastic 
theory adopted, discloses a deep and fascinating feature 
of curved crystals by requiring isometric surfaces (as the 
catenoid and the helicoid or the Scherk surfaces) to sup- 
port the same crystalline order. From the other hand this 
observation can be used to set a bound in the complete- 
ness of our theory and the accuracy of linear elasticity in 
describing the physics of defects in curved crystals. We 
hope more work will be devoted to the understanding of 
these phenomena in the future. 



Acknowledgments 

This work was supported by the NFS through Grant 
No. DMR-0219292 (ITR) and through founds provided 
by Syracuse University. 

APPENDIX A: EVALUATION OF THE 
FUNCTIONS G L (x) AND T(x) 

In this appendix we discuss the calculation of T(x) in 
Eq.(17). Consider the solution of the Green equation 

A 2 G 2L (x,y) =6(x,y) x,yeF (Al) 

with homogeneous boundary conditions. In integral form 
this solution can be written: 

G 2L (x,y) = J d 2 zG L (x,z)[G L (z,y) - H(z,y)] (A2) 

where Gi(a;,y) is the solution of the Green-Dirichlct 
problem Eq.(12) and H(x,y) is the reproducing kernel 
of Eq.(23). As noted in Sec. II A, one can map the 
paraboloid P onto the unit disk of the complex plane 
and then employ the appropriate planar techniques (i.e. 
image charges). In general any simply connected two- 
dimensional Ricmannian manifold with a C°°— smooth 
metric ds can be equipped with a set of local isothermal 
(or conformal) coordinates (cc, y) such that the metric is 
represented in the form ds 2 = w(x, y)(dx 2 +dy 2 ) for some 
positive conformal weight w. In two dimensions this local 
system of isothermal coordinates can serve as a conformal 
chart for the unit disk D on the complex plane. Calling 
z = ge 1 ^ the new metric will be 

ds 2 =w(z)(dg 2 + g 2 dcj) 2 ). (A3) 

The conformal factor w{z) can be found by equating the 
metric (A3) with the original one. At this point it is 
worth treating the problem in a slightly more general 
form. Consider the case of a generic surface of revolution 
of the form: 

x = £(r) cos 4> 

y = £(r) cos 4> , (A4) 
z = r)(r) 

with r € [0, oo [ and <fi £ [0, 2ir] . The metric of the surface 
(A4) will be: 

ds 2 = Edr 2 + 2Fdrd<j) + Gd<j) 2 , (A5) 

where 

E = C 2 + rf 2 

F = (A6) 
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are the coefficients of the first fundamental form of the 
surface (A4). Equating (A5) and (A3) one finds imme- 
diately: 



>(Q) 



n 2 



with q and r related by the differential equation: 
dg.fE 

whose solution is given by: 



g = exp(r J dr^fEjG\ 



(A7) 



(A8) 



(A9) 



The sign of the exponent and the integration constant 
in (A9), can be tuned to obtain the desired scale and 
direction of the conformal map. It is easy to show that, 
in the new coordinates, the Laplace operator takes the 
form: 

A = w~ 1 A z , (A10) 
in which A z is the Laplacian in the Euclidean metric: 

(All) 



7 : 



1 

Q 2 



With this new set up, the Green equation (Al) can be 
easily written in the form: 



A., 



where £ = r 1 e l 



L A 2 G 2i (z,C)=<5(z,C), 



(A12) 



is a second generic point of the com- 
plex plane and S(z, Q has the same meaning as in Sec. 
II A with respect the Euclidean metric 7. The differen- 
tial operator Aw -1 A is known in analysis as weighted 
biharmonic operator. Analogously, the Green-Dirichlct 
problem (12) becomes: 



so that: 



AGl(z, •) = 5(z, •) zeD 
G L (z,-)=0 zedl 



G L (z,C) = — log 



1-zC 



(A13) 



(A14) 



It must be stressed that so far we didn't explicitly use 
the geometry of the paraboloid. What has been said, 
therefore, holds for any surface of revolution which can 
be conformally mapped onto the unit disk. Furthermore, 
as anticipated in Sec. II A, the conformal distance g, 
which completely embodies the geometry of the surface, 
depends only on the coefficients E and G of the first 
fundamental form. 

In the particular case of the paraboloid we have: 



A 



re 



1 + vT + K 2 r 2 ' 



(Al5a) 



A = 



1 + VI + k 2 R 2 



i?cxp( v / l + K 2 i? 2 ) 



(A15b) 



To obtain the expression for T(x) given in Eq.(17) we are 
left with the task of calculating the integral: 



r s (r) = r S)1 (x)-r s , 2 (r), 



(A16) 



where 

r a ,i(aO = ^ / dcj>'dr'^gK(r')log\z-(\, (Al7a) 
Vs ^ x) = h (<WdT> y /gK(r')log\l-z(\. (A17b) 
For this purpose we can use the expansion: 

00 1 / \ n 

log|z-C| =loge> - V- — cos n8<P, (A18) 

where q > (f?<) represents the largest (smallest) modulus 
between z and £, while S(f> = <f> — <$>' . The factorization 
of the angular variables in Eq.(Al8), together with the 
pure radial dependence of the Gaussian curvature and 
y/g, makes the angular dependence of T s i vanish, so that 
we have: 



r. ll (r)=lo ge (r) f dr> ^ 
Jo (1 + K z r u 



dr' 



(1 + K 2 r' 2 )2 
which integrated by parts gives: 



! )3 



loge(r'), (A19) 



r ail (r)=log 



ae 



1 + Vl + n 2 r 2 



(A20) 



Using an expansion similar to (A18) it is also possible to 
prove that 

00 

log |1 - zC| = - V - (BQ') n cos^, (A21) 

z — J m 



which integrated over the surface of the paraboloid gives 
F s 2 = 0. This last conclusion, combined with Eq. (A20), 
yields Eq.(18). 



APPENDIX B: THE HARMONIC KERNEL 

In this appendix we give the exact expression for the 
harmonic kernel appearing in Eq.(23). This expression 
has been found by Shimorin in the more general case of 
the calculation of the Green function for the weighted 
biharmonic operator with radial weig ht w(z) = w(\z\ 2 ) 
[16]. As in the previous appendix we call z — p(r)e 1 ^ 
and ( = p'(r')e 1 ^ two points of the complex plane that 
arc images of the points (r,<j>) and (?',(/>') of P under 
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FIG. 8: CPU- Time and speedup (i.e. the time employed by 
n processors to accomplished a given number of iterations 
divided by the time employed by a single machine to achieve 
the same task). 



the conformal transformation (16). The harmonic kernel 
H(z, can be written in integral form as: 

"o-o—CSf.**"^*)- (B1) 



in which: 



(B2) 



n>0 n<0 

where the coefficients c n are given by: 



Cn = 2 [ ds^ Q 2n {s). (B3) 
Jo 



APPENDIX C: OPTIMIZATION VIA PARALLEL 
DIFFERENTIAL EVOLUTION 

Many of the techniques proposed to determine the 
crystalline structure of systems of interacting particles, 



as in the Thomson problem, are based on local opti- 
mization procedures such as steepest descent, conjugate 
gradient and the quasi-Ncwton method. Such methods 
belong to the class of line-search algorithms for multidi- 
mensional non-linear programming problems. They can 
be described, in general, as a sequence of line minimiza- 
tions of an objective function along a set of directions 
that are generated differently in different algorithms. Be- 
sides the well known local convergence properties of these 
methods, they are generally unable to locate the global 
minimum since they inherently approach the closest local 
minimum for a given set of initial conditions. 

To avoid the misconvergence problem described we 
adopt the Differential Evolution (DE) algorithm of Storn 
and Price [28]. This algorithm, which has been success- 
fully applied to several optimization problems in engi- 
neering [29] , belongs to the family of evolutionary algo- 
rithms which are considerably faster than other stochas- 
tic optimization methods, such as simulated annealing 
and genetic algorithms, and more likely to find the cor- 
rect global minimum. These methods heuristically mimic 
biological evolution by implementing natural selection 
and the principle of "survival of the fittest" . An adaptive 
search procedure based on a population of candidate so- 
lutions is used. Iterations involve a competitive selection 
that drops the poorer solutions. The remaining pool of 
candidates are perturbed (or mutated) in order to gener- 
ate trial individuals and then rccombincd with other so- 
lutions by a swap of the components. The recombination 
and mutation moves are applied sequentially; their aim is 
to generate new solutions that are biased towards subsets 
of the search space in which good, although not necessar- 
ily globally optimized, solutions have already been found. 

An essential feature of Differential Evolution is the es- 
tablishment of genetic diversity, which helps to maximize 
the probability of finding the true global minimum and 
to avoid misconvergence. One begins with a large popu- 
lation of individuals uniformly distributed in the search 
space. A good choice, in practice, is to choose the number 
of individuals to be an order of magnitude more than the 
number of variables in the problem. The price one pays 
is a dramatic slowing down of the algorithm when ap- 
plied to large scale optimization. Considerable effort has 
therefore been made in the past ten years to develop par- 
allel implementations of evolutionary algorithms aimed 
at reducing the overall time to completion of the task 
by distributing the work on different processors working 
in parallel. More recently some researchers have con- 
jectured that some parallelizations of a task improve the 
quality of the solution obtained for a given overall amount 
of work (e.g. emergent computation). 

The Island Model is a popular choice among paral- 
lelization strategies and is implemented within a message 
passing model. It consists of dividing the initial popu- 
lation into several sub-populations and letting each of 
them evolve independently on a single machine for a pre- 
determined number of iterations (called the epoch) . The 
exchange of genetic information is promoted by swapping 
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individuals between different sub-populations at the end 
of each epoch. In the present work the migration strat- 
egy consists in swapping the best individual of each sub- 
population with a randomly selected individual on an- 



other island with the ring topology chosen for the connec- 
tivity between islands. This choice allowed us to achieve 
a substantial reduction of the CPU time and a linear 
speedup (see Fig. 8). 
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